约束聚类-多元回归树及重要判别变量识别及R操作

前文已简介了决策树，其基本思想是对预测变量进行二元分离，从而构造一棵可用于预测新样本单元所属类别的树。常规的决策树中，响应变量只包含单变量的分类或定量变量，因此也称为一元回归树。

多元回归树（multivariate regression tree，MRT）是一元回归树的拓展，可用于描述和预测两个变量矩阵间的关系，其响应变量是多变量结构的（De’ath, 2002）。MRT是多元回归的一种形式，响应变量可以通过解释变量进行解释，并且可以通过解释变量进行预测，是一种约束聚类方法。MRT和冗余分析（RDA）存在相似之处，即都通过解释变量对响应变量实现“约束”，二者相比MRT重在关注解释变量的预测能力大小，RDA重在关注解释程度的大小。

多元回归树概述

MRT由两步程序组成，数据的约束划分和分类结果的交叉验证。

其计算过程可概括如下。

（1）随机将对象分为k组，选出1组作为验证组，其余k-1组重新混合并通过约束划分（以多元回归为基础，通过解释变量拟合响应变量）建立回归树，通过最小化组内平方和（最小化类内对象的差异）识别分类；

相对误差（relative error，RE）是回归树中所有分支的组内平方和与响应变量的总体平方和的比值，即响应变量总方差不能被回归树解释的比例，RE等同于1-R²，R²则代表了回归树的解释能力程度；

（2）重复（1），即依次剔除1组数据，共重复（1）k-1次；

（3）最终获得k个回归树，将验证组的对象分配到每个分类水平中，并由下式计算交叉验证相对误差（cross-validation relative error，CVRE）：

y_ij(k)为验证组k中的每个观测值，ŷ_ij(k)为验证组k中每个观测值的预测追（分类组形心），分母代表响应变量的总体平方和；因此CVRE也定义为验证组未能被树解释的方差与响应变量总方差的比值，值越低代表模型预测性能越好；

（4）裁剪回归树：保留具有最小CVRE或保留CVRE在最小(CVRE+1)个标准差范围内的最小分类水平的回归树；

（5）重复n次（1）-（4），获得误差估计值；

（6）保留具有最小CVRE或保留CVRE在最小(CVRE+1)个标准差范围内的最小分类水平的回归树。

在实际的回归树选择中，最好综合考虑RE、CVRE以及树规模的大小，尽可能选择RE、CVRE以及树规模都低（预测精度高且简约）的树模型。

MRT最早应用在生态学研究中。MRT通过对物种多度数据的重复分割形成聚类簇，每个分割都由基于环境值的简单规则定义，通过选择物种差异的度量，可将物种组成的任何方面与环境数据相关联。聚类簇及其对环境的依赖关系通过树状图表示，每个簇代表了物种集合，其环境值定义了与其相关的栖息环境特征。MRT可用于分析复杂的生态数据，包括失衡、带缺失值、变量之间的非线性关系和高阶相互作用等。此外，MRT还可用于获得指示环境的关键物种，以及预测仅获得环境数据的地点的物种组成。

接下来就展示MRT在R中的实现方法，以及如何获得重要的变量。

R包mvpart的多元回归树

在R中，MRT可使用mvpart包执行。

mvpart包无法使用install.packages()或bioconductor直接安装，因此可能需要折腾一下。

##安装
#参考自
# https://stackoverflow.com/questions/29656320/r-mvpart-package-any-option-to-use-in-r-3-1-x#opennewwindow

#首先安装 Rtools，点击以下链接
# https://cran.r-project.org/bin/windows/Rtools/

#Rtools 安装好之后，再安装 mvpart 等包
install.packages('devtools')
devtools::install_github('cran/mvpart') #计算多元回归树
devtools::install_github('cran/MVPARTwrap') #辅助分析

#要是在已经安装了 Rtools 之后，再安装 mvpart 等包时还提示“ had non-zero exit status”之类的
#重新再安装一个新的 R 后，再安装 mvpart 等包

#加载
library(mvpart)
library(MVPARTwrap)

数据集

数据集“spider”，包含28项观测值，其中前12列为不同蜘蛛物种的丰度观测值，后6列为环境变量的测量值。所有数值已经过某种标准化处理，并非原始测量值。

#数据集，详情 ?spider
data(spider)
head(spider)

接下来使用该数据集展示R包mvpart的MRT方法。

注：MRT中，平方和的计算具有欧式空间属性，某些响应变量的原始值可能不适合直接使用。因此使用自己的数据分析时，需要考虑数据本身的特点，原始数值是否需要进行某种转化。例如，量纲不同的变量通常需要标准化为均值0、标准差1的结构；原始物种多度数据通常需要预转化（如弦转化）才可应用欧式空间等。

构建MRT

生成一个被解释变量（环境组成）约束的响应变量（物种多度）的多元回归树。

#数据集，详情 ?spider
data(spider)
head(spider)

#构建 MRT，详情 ?mvpart
#~. 意为使用所有环境变量，等同于 ~herbs+reft+moss+sand+twigs+water
species <- spider[ ,1:12]
env <- spider[ ,13:18]

set.seed(123)
spider_mvpart <- mvpart(data.matrix(species)~., env, xval = nrow(species), xvmult = 100, xv = 'pick', pca = TRUE)
spider_mvpart

该图为相对误差（relative error，RE）和交叉验证相对误差（cross-validation relative error，CVRE）与分类组数（size of tree，即树的分支数）的关系图。误差图为交互式的，允许自定义选择合适的分类方案，点击误差图的特定位置即可展示回归树结构。

一般来讲，RE或CVRE会随着树的分支数的增加而减少，开始时下降比较明显，但到某范围后将变得平缓，甚至可能会有所上升。图中红点位置表示具有最小CVRE的分类方案，即代表了具有最好预测能力的树，上方绿色条形代表了交叉验证迭代的次数。

但实际中，并非一定要选择最有最小的误差值，还要避免树过于复杂。我们看到，尽管最小CVRE在树的分支为5时出现，但在3或4时即达到了一个相对较低的平缓值，其后随着分支的增加误差下降并不明显。由于数据集本身并不大，简约起见不妨就选择分3组的回归树，它的预测能力与具有最小CVRE的回归树接近（精度相似），但更便于观测。

回归树中展示由解释变量预测的群落结构，每个判别节点处为物种响应的环境梯度，基部终端节点处的不同颜色的条形图代表不同物种，高度代表物种的丰度。数值n代表了分类至该组中的样方数量。

回归树的底部还展示了一些常规的统计信息，误差（error）、交叉验证误差（CV error）、标准误差（SE）等。

继续点击图中的任意节点，可展示该数据集物种丰度组成的主成分分析（PCA）结果。

PCA图中的点表示样方，并按分类颜色标注出。12个向量代表了12个物种，延伸方向即为这些物种丰度增加的方向。详见PCA排序图说明。

MRT更多细节概要

通过summary()，可查看更多细节，如各个节点的特征等。

summary(spider_mvpart)

例如CP表中汇总了不同大小的树对应的预测误差，可用于辅助设定最终的树的大小。

CP表中各部分的解释，可参考前文决策树。

#提取需要的结果用于进一步观测，例如获取 CP 表
#names(spider_mvpart)
spider_mvpart$cptable

获取回归树对样方的分类详情，上文中已手动选择聚为3类。

#例如查看样方分组
spider_class <- spider_mvpart$where
names(spider_class) <- rownames(spider)
spider_class

获得样方划分分组后，可进一步作图展示各组中的物种丰度信息，如堆叠柱形图、饼图、排序图等，区分组间差异，不再多说。

R包MVPARTwrap的辅助功能

MVPARTwrap包提供了展示回归树的更多细节的功能。

MRT()主要结果包括：回归树中的成员分配，每个分类节点的判别物种，节点的R²（1-相对误差SE），物种解释方差贡献，回归树各支解释方差等。

#MVPARTwrap 包的辅助功能，详情 ?MRT
#该包的安装方法同 mvpart，已在上文提及
library(MVPARTwrap)

spider_mvpart_wrap <- MRT(spider_mvpart, percent = 10, species = colnames(species))
spider_mvpart_wrap

summary(spider_mvpart_wrap)
plot(spider_mvpart_wrap)

总结项比较多，以下仅选择部分结果用于展示。

每个物种在每个节点对方差解释的贡献度，R²（1-相对误差SE）越高代表了对该节点的划分越重要，即暗示该物种可作为该节点的判别物种来看待。

指示种（indval species）显示了每个节点的判别物种，即对方差解释贡献最大的物种。指示值越高表明该物种越重要，显著性由置换检验获得。

判别的样方所属的聚类群，上文中选择聚为3类。

参考资料

DanielBorcard, FranoisGillet, PierreLegendre, et al. 数量生态学:R语言的应用（赖江山 译）. 高等教育出版社, 2014.

De’ath G. Multivariate regression trees: a new technique for modeling species–environment relationships. Ecology, 2002, 83(4):1105-1117.

友情链接